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We study quantized solutions of the Wheeler de Witt (WdW) equation describing a closed 
Friedmann-Robertson-Walker universe with a A term and a set of massless scalar fields. We show 
that when A <C 1 in the natural units and the standard in-vacuum state is considered, either wave- 

function of the universe, T, or its derivative with respect to the scale factor, a, behave as random 
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quasi-classical fields at sufficiently large values of a. The former case is realised when 1 -C a -C esK, 

2 

while the latter is valid when a ess. The statistical r.rn.s value of the wavefunction is pro¬ 
portional to the Hartle-Hawking wavefunction for a closed universe with a A term. Alternatively, 
the behaviour of our system at large values of a can be described in terms of a density matrix 
corresponding to a mixed state, which is directly determined by statistical properties of T. Similar 
to T, the density matrix can be considered as c-number valued in the position and momentum 
representations. The probability distribution to find a universe with particular values of the scale 
factor and field amplitudes following from this density matrix is again proportional to that of the 
Hartle-Hawking wavefunction, while the probability distribution over field velocities is non-trivial 
and different from what follows from the Hartle and Hawking formalism. 

We suppose that the same behaviour of \F can be found in all models exhibiting copious production 
of excitations with respect to the out-vacuum state associated with classical trajectories at large 
values of a. Thus, the third quantization procedure may provide a ’boundary condition’ for classical 
solutions of the WdW equation. Contrary to the previous proposals, in our case two equivalent 
descriptions of this classical solutions are possible. Either <F can be regarded as a stochastic classical 
quantity or the system can be viewed as being in a mixed state defined over classical solutions to 
WdW equation. 

Keywords: Quantum Cosmology, Third Quantization, closed FRW Universe, Lambda term 


I. INTRODUCTION 

One of the possible approaches to the problem of quantum gravity is based on the canonical quantization of 
Einstein equations and analysis of properties and solutions of the emerging Wheeler de Witt (WdW) equation for 
the wave functional of the Universe, 4* [1], see e.g. [2] for a review and further references. When *F is considered 
as a c-number value some long standing issues arise. First, a direct probabilistic interpretation is hampered, second, 
there is an ambiguity related to initial (or, boundary) conditions for solutions to WdW equation. Over the years 
several proposals for a possible choice of initial conditions have been made, notably the ones leading to the tunneling 
wavefunction of Vilenkin [3] and ’no-boundary’ wavefunction of Hartle and Hawking [4]. 

Some problems are alleviated in an approach, where itself is treated as operator valued, IF —>■ >F. This is based 
on observation that the WdW equation has a formal structure of a hyperbolic equation with a variable determined by 
the volume element of spatial hypersurfaces playing a role of a ’time’, which can again be canonically quantized. This 
procedure is called ’the third quantization formalism’, see e.g. [5]. It was explicitly realised in minisuperspace models, 
where the metric was fixed by the condition that it is of Friedmann-Robertson-Walker type, and matter degrees of 
freedom were modeled as a set of scalar fields, see e.g. [7], [6]. 
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In the minisuperspace approach the wave functional becomes a function depending on the scale factor, a, and field 
amplitudes, and WdW equation is reduced to a Klein-Gordon equation with second derivatives over the scale factor 
and the amplitudes entering it with opposite signs and, a potential being, in general, a function of these variables. 
When the logarithm of the scale factor is chosen as a time variable the differential part of WdW equation takes 
the standard form of the d’Alembert operator in a certain factor ordering scheme assumed from now on. Close to 
singularity this variable tends to minus infinity and it may be shown that the potential tends to zero. In this limit the 
quantized WdW equation is the familiar Klein-Gordon equation for a massless quantum scalar field, and, therefore, 
the standard vacuum state for this field can be chosen as a quantum state of the system, e.g. [6]. Hereafter, we 
call it the in-state having in mind the analogy with formalism of quantum fields in curved spacetimes, see [8]. On 
the other hand, in a class of models, where classical dynamics with large values of scale factors is possible, there is 
another natural vacuum state associated with a set of quasi-classical solutions of WdW equation called hereafter the 
owt-state. In general, this state does not coincide with the in-state, thus the latter contains excitations with respect 
to the former one. This effect is interpreted as ’creation of universes from nothing’ in the framework of the third 
quantization procedure [30]. Technically, this means that a positive frequency solution defined with respect to the 
in-state is a mixture of positive and negative. Care should be taken with this terminology. Rather, as in the general 
case of creation of particles by an external field, see e.g. [8], we are talking about amplification of vacuum fluctuations 
of the quantized wavefunction by evolving with the scale factor gravitational field.frequency solutions corresponding to 
the out-state, with the Bogolyubov coefficient j3 in the decomposition of positive frequency solutions over the negative 
frequency ones being nonzero. The amount of produced universes with given constants of motion corresponding to a 
particular quasi-classical out mode is proportional to the square of absolute value of the /3, summed over all modes 
defined with respect to the ?'n-state. 

In this Paper we would like to propose another interpretation of the emergence of ’classical’ properties of the 
system on hand in the framework of the third quantization formalism. Namely, we show that in a range of sufficiently 
large scale factors the quantized wavefunction can behave as a random classical variable. We consider a simple WdW 
minisuperspace model of a closed FRW universe with a A term and n homogeneous massless scalar fields, pi, 2 < * < n. 
It was shown that when A <C 1 in the natural Planck units, typical values of /3 are exponentially large, j3 oc e* [7]. 
On the other hand, the WdW equation of this model is quite similar to those used to describe the dynamics of test 
quantum fields in a natural in -vacuum state in inflationary models, see e.g. [II], [9], [10]. As is known the latter 
problem also exhibits a copious ’particle’ production with respect to a suitably chosen oid-state. Also, a large scale 
part of these fields essentially behaves as a random classical field, e.g. [11], [13], [12]. In this case the creation and 
annihilation operators in the decomposition of the test field over the normal modes can be treated as classical random 
variables obeying Gaussian statistics. There are several possible tests of this (quasi)classical behaviour, in particular 
based on investigation of the evolution of the in-vacuum state in the Schrodinger representation (e.g. [14], [13]), or the 
use of the Wigner function, e.g. [15], [13]. In particular, the in-vacuum should be strongly squeezed with dispersion 
of a canonical variable manifesting a quasi-classical behaviour being much larger than a typical one following from 
the uncertainty principle (e.g. [14]). The Wigner function in the regime of quasi-classical dynamics is approximately 
reduced to a phase density distribution of a bunch of classical trajectories. 

We show that both these criteria are satisfied in our model in a certain average sense, for modes giving the main 
contribution to the expectation values of interest [31]. Therefore, these expectation values can be treated as statistical 
averages of classical quantities. However, an important difference with the test field case consists in the fact that in our 
case mode amplitudes, and, accordingly, T, are approximately quasi-classical only when l<a< eax. When a esx 
the mode momenta are quasi-classical, and, therefore, it is the derivative of T over a, which exhibits quasi-classical 
behaviour. Thus, the third quantization procedure together with a natural choice of the m-state may be used to define 
an initial condition for a classical wavefunction ’F. Unlike the known proposals for the initial conditions in our case 
d' is essentially a random quantity. Interestingly, the rms value of our wavefunction, \/< Tfl/* >, is approximately 
proportional to the Hartle and Hawking expression in our model, although it is not clear to us, whether this fact is 
due to a coincidence, or it may be of a generic nature. 

Alternatively, when 1 C a < ea we can use a density matrix with c-numbered matrix elements instead of the 
classical wavefunction. Thus, in our model in this regime the quantum state of the Universe is a mixed one. Its 
diagonal elements in the position representation are equal to < ’F’F* >, while it has a non-trivial structure in the 
momentum representation giving a probability to find universes with different values of pi. 

The structure of the Paper is as follows. In Section II we introduce basic definitions and equations. In Section III 
we obtain asymptotic solutions to the WdW equation in the limit A —> 0 based on a procedure involving the Wentzel- 
Kramers-Brillouin-Jeffreys (WKBJ) technique and the use of different solutions with common ranges of validity. In 
Section IV we discuss values of |/3| 2 and the behaviour of the Wigner function and the vacuum state in our model, 
showing that they all indicate quasi-classical dynamics of the variables of interest in the sense explained above. In 
Section V an explicit form of ’F valid for values of a 1 is obtained and its averaged value is calculated. In the same 
Section we derive expressions for the density matrix in the momentum and position representations and discuss its 
properties. Finally, we present our conclusions and discuss our results in Section VI. 

We use below the natural Planck system of units setting Planck and gravitational constants as well as the speed of 
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light to unity. 


II. BASIC DEFINITIONS AND EQUATIONS 

In what follows we are going to consider the quantum dynamics of a FRW universe with positive spatial curvature, 
having a cosmological term A, and n massless scalar fields ip t , where i = 1, ..n. In this case the WdW equation takes 
the form (e.g. [6], [18]) 


0 2 A 

^2 - A„ + - exp (6f) - exp (4 1) ) T = 0, 


( 1 ) 


where t = In a and a is the scale factor [32], A„ = 5 ZT=i and is the wavefunction. In agreement with the third 
quantization procedure we treat the wavefunction as an operator obeying the standard commutation relations: 


r), = iS n (x - x), 


( 2 ) 


where [..] is the field commutator, x represents n-dimensional vector with components ifi, 5 n (x) is n-dimensional 
Dirac delta function and the dagger stands hereafter for Hernritian adjoint. 

Equation (1) is formally a Klein-Gordon equation for a free quantum field having a time dependant potential 


V(t) = ^ exp (6f) - exp (4f), 

and the associated Lagrangian and Hamiltonian can be written in the form 


L = J d n x (^ - VM^ , H = J d n x (PP* + , + F'F’F t ) 


( 3 ) 


( 4 ) 


where commas stand for partial differentiation over tpi, summation over repeating Latin indices is implied from now 
on, and P, P^ are the canonical momenta. From the Hamilton equations we easily obtain 


P = iU t 
m ’ 


=P 


( 5 ) 


Note that the variable t plays the role of a time-like variable in equation (4), and, therefore, it will be referred to 
as ’time’ below. In order to avoid confusion, let us point out, however, that the proper time is clearly absent in the 
WdW equation due to its invariance with respect to a choice of the lapse function. 

Solutions to (1) can be represented in the standard form 


>F = 


J d n k ( 


U u e ik - x d k + U:e- ik - x bl, , 


( 6 ) 


where k is an n dimensional vector, • stands for the scalar product, * is the complex conjugate and u> = y/k ■ k. The 
mode functions U, u satisfy the equation 


U u + (w 2 + V)U U = 0, 


( 7 ) 


where dots stand for differentiation over the tinre-like variable t. Note that we clearly have U u = U- 
Assuming that the mode functions are normalised according to the condition 


TI TJ* - U U* — _ 

U “ U U1 — ( 27r )n ’ 


( 8 ) 


the operators dk and bj- obey the standard commutation relations for the creation and annihilation operators 

[d k ,d\,)=5 n {k-k), [b k , &[.,] = 5 n {k - k) (9) 

with other commutators being equal to zero. 

It is worth noting that near the singularity when t —> — oo the potential V tends to zero and equation (1) formally 
describes a massless free scalar field in an effective Minkowski spacetime. Thus, when the field mode U u is an 
eigenfunction of the timelike Killing vector in this spacetime, and, accordingly, oc e~ lult , the vacuum state |0) 
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defined in such a way that dk |0) = 0 and bk |0) =0 for all k is the standard vacuum state for a massless scalar field 
in this asymptotic limit. We use this state as the field state in our analysis below. 

For our purposes it is sometimes convenient to use another representation for the field W through a Fourier transform 

* = J d n ke ikx ^ ul , *„ = -^= (U u (ci, k + *C 2 ,fc) + U*{c\ k + *4 >fc )) , (10) 

where we introduce new creation and annihilation operators 


Cl,k — ^(h/c + fc)? 




( 11 ) 


It is easy to see that these operators obey the standard commutations. Also, it is evident that the vacuum state |0) 
defined above is also a vacuum with respect to these operators. 

The Hamiltonian (4) can be expressed in terms of 46% as 

H = (2tt)" J d n k (P U P' + (cu 2 + TOM!) , P u = JU*. (12) 


It can be brought to a standard form by separating into the real and imaginary parts as 41 w = . 1 (qi u + iq 2 w), 

v' 2 ( 27r ) n 

where q± tU and q 2 lU are Hermitian and commute with each other. We have 



+ (w 2 + H)g 2 >w ) , 


Pa,co 


d_„ 

di Qa ’“ 


(13) 


The expression (13) tells that the problem can be formulated in terms of an infinite set of oscillators with a time 
dependent frequency. Based on the analogy with the oscillators we introduce other time-dependent “creation and 
annihilation operators”, dk (t) and d\. (t ), according to the rule 


4(1)^^ + ^=^, (14) 

where 

fi(to) = Pu; 2 + V(to), (15) 

the potential V is assumed to be taken at a particular fixed moment of time, to, and, therefore, f2(to) does not depend 
on time t. Also, for simplicity we assume hereafter < 0 is such that w 2 + V(t-o) > 0, and, accordingly, fl(t 0 ) is real. By 
construction the operators (14) obey the standard commuting relations. Thus, they can be related to the operators 
Cfc by a Bogolyubov transformation 

dk(t) = a u Ck + P u c\, |c^| 2 - |/3 W | 2 = 1- (16) 

The explicit form of the Bogolyubov coefficients follows from expressions (10-14): 

= + = + (17 > 

The normalisation condition (8) tells that the Bogolyubov coefficients do obey the second equality in (16). 

The vacuum state defined by the condition dk(to) |0) ad = 0 defines so-called “adiabatic” vacuum state |0) ad (see e.g 
[8]). Note that, by definition, this state does not depend on time. 

Of special importance are two particular values of to, to ±oo. 

In the former case the adiabatic vacuum state provides a natural out-vacuum state for our problem. Thus, when 
fik(to ~> oo) ^ 0 the state |0) contains excitations above the out-vacuum usually interpreted as “creation of particles” 
(or universes, in our case) from ’nothing’. 

The latter case corresponds to setting V(to) in (15) to zero. Clearly, in this situation relations (16) and (17) 
simply determine the evolution of Ck and c\ in the Heisenberg representation. On the other hand, in the Sclirodinger 
representation these relations can be used to determine the evolution of the wave functional, see e.g. [13] and Section 
IV B below. 
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III. AN ASYMPTOTIC ANALYTIC SOLUTION OF THE WDW EQUATION 

In general, solutions to equation (7) can be obtained analytically only when either the first or the second term 
in the expression for potential V in (3) is discarded. This corresponds to neglecting the influence of the A term or 
spatial curvature, respectively. When both terms in equation (3) are retained it can either be solved numerically 
or an approximate solution can be looked for using appropriate asymptotic methods. Here we consider the second 
possibility in detail in the limit A 1, which may be appropriate for inflationary models, where the value of the A 
term is generally much smaller than its natural Planck value A = 1. Numerical solutions will be used to validate our 
analytic approach. 

Additionally, for simplicity we are going to consider sufficiently small values of the frequency cu, namely, in our 
analytical work we formally assume hereafter that to -C uj cr it, where LU C rit = y/—V m in{tmin) = and Fmin(tmm) = 

— , tmin = t; In are the minimal (negative) value of the potential V and the corresponding value of the variable 

t , respectively. For such values of u we can use the WKBJ approximation, assuming that a solution to (7) is 
approximately proportional to 

-4= e iS{t) , S = ± 

both in the classically forbidden region corresponding to an intermediate range of t, where w <C — V and in the 
classically allowed regions corresponding to small values of t, where V —> 0, and large values of t , where V ui. In 
the forbidden region the phase S is purely imaginary, and therefore a general solution to (7) has the form of a sum 
of growing and decaying exponents multiplied by the pre-exponential factor, while in the allowed region S is real and 
the solution to (7) has oscillatory behaviour. It is very important to note that the WKBJ approximation breaks down 
when the time t is close to = i In such that U(t*) = 0. However, in the vicinity of this moment of time one can 
simplify (7), find an appropriate exact solution of the simplified equation and match it to the WKBJ solution. 


U)‘ 


+ V, |S|»1 


(18) 



FIG. 1: We show the potential V(t) together with four overlapping regions, where various approximations discussed in the text 
are possible. The potential is calculated for A = 3 * 10 —3 , square of the mode frequency is taken to be co 2 = 0.5. Since the 
precise positions of boundaries of these regions are ambiguous they are shown schematically. Vertical solid, dashed, dotted, 
dot-dashed lines show the positions of boundaries of regions 1), 2), 3) and 4), respectively. 


Thus, the whole time interval —oo < t < +oo can be subdivided into four overlapping regions: 1) the region of 
sufficiently small t, where the curvature term dominate over the term proportional to A in the expression for the 
potential, 2) the classically forbidden region, 3) the region close to the moment t = i* and 4) the classically allowed 
region. These regions together with the potential V(t) are schematically shown in Fig. 1. Matching solutions in 
all these regions and using the normalisation condition (8) we can find a solution for U,^ approximately valid when 
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ui <C uj C rit • As we shall see below, a comparison with numerical results shows that this solution is qualitatively valid 
even when oj < LO cr i t . On the other hand, an approximate solution in the opposite limit ui to cr it can be obtained 
by setting the term determined by the curvature in (3) to zero. Its form is well known, see e.g. [2] and unimportant 
for us. 

Let us consider the behaviour of U u by turns, starting from region 1). 

In this region the explicit form of the solution can be written as 


U u 


CLJ_ i4f 



r(i-f) 

) n uj 


(19) 


where I v (x) is the modified Bessel function of the first kind, the coefficient C u is determined by the normalisation 
condition ( 8 ) and T(a:) is the gamma function, see also e.g. [17]. It is easy to see that when e 2t <C 1 from equation 


(19) it follows that U u 




=e lut , which is just the standard positive frequency solution in (n + l)-dimensional 


Minkowski spacetime normalised according to ( 8 ), see e.g. [18]. 

Below we shall need an asymptotic form of the modified Bessel function of purely imaginary index at large values 
of its argument. As discussed in e.g. [19] to this purpose it is convenient to consider the modified Bessel function of 
second kind K v (x) and an additional function L v {x) defined according to the rule 


7r 77T 

K A X ) = '■ 7 —T ( T -Ax) - h{x)) , L v {x) = ———T (I-v(x) + I v {x)) . (20) 

Note that both functions are real when their argument is real and the index is purely imaginary. As shown in [19] 
when x —> 00, Li ll (x) « while Ki /J/ (x) has the well known asymptotic form: K^x) « \f^e~ x in the 

same limit. Using these expressions and considering sufficiently large values of t we obtain from (19) and (20) 



Note that although the last term in the brackets is exponentially small at large t it is needed to be retained for the 
asymptotic solution ( 21 ) to satisfy the normalisation condition ( 8 ). 

The solution (21) can be matched to a WKBJ solution of the form (18) in a time interval within region 2), where 
on one hand the time t is small enough such that t < t m i n , and on the other hand, it is sufficiently large for the 
condition to <C yJ—V to be fulfilled. 

Before doing so let us discuss the WKBJ phase, S, in (18). In general, it can be represented in the form 


S = ±i j ^ \/P(x), P(x) = 4x 2 - ^Aa ; 3 - 


CJ 


( 22 ) 


where we introduce a new independent variable x = Since the polynomial P(x) is of the third order explicit 
integration over x in (22) is rather cumbersome. It becomes trivial, however, when we set ui = 0. In this case 
P(x) (x x 2 , and the integration gives S = y^(l — 2A£) 3 / 2 . In order to avoid unnecessary complications let us take 
into account that in the limit u> <C LO cr it the term w 2 in (22) may play a role only when t ~ i*. To take this fact into 
account we can consider a modified polynomial, P m od , in (22) fixed by the following conditions: it is of third order, 
Pmod, oc x 2 , Pmod, = P when w = 0 and P mo d(t*) = 0 . By doing so we obtain P mo d = x 2 {4 — f A 2 w 2 — |Ax) and 

S = q=x(l - ^if 1 - ^) 3 / 2 - 

Taking this into consideration it is easy to see that the matching procedure in region 2) gives in the leading order 


EA 


CL 


A / 4 


(ei 




a/ 2 '. 


. . . /7TW 

i sinh e 




3/2> 


(23) 


where 


<j>- = 1 — 


A 2 w 2 

9 



(24) 


2 a 2 

Now let us consider region 3), where t ~ i*. For that we assume that relative difference (f* — -t)/t* is small 

and simplify equation (7) retaining only leading order terms with respect to the relative difference. It is convenient 

to introduce a natural variable z = (t* — - t) to see that after the simplification equation (7) is reduced to 

the Airy equation 



zU k = 0 . 


(25) 



Solutions to (25) are matched to (23) in the region, where both conditions (t 
fulfilled. This gives the solution in the form 
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-t)/t* <C 1 and z 1 are 


1 

U u = C^ 6 (jj (2e*Ai(z) + isinh e~*Bi(z)^ , (26) 

where Ai(z) and Bi(z) are Airy functions of the first and the second kind, respectively. 

The solution can be analytically continued to the region z < 0. Using a standard result and assuming that 
z\ = — z 'A> 1 we obtain 



Finally, we can match solution (27) to the solution of the form (18) in the classically allowed region 4) using 
the same technique as the one leading to expressions (26) and (27) and taking into account that in this region 
S = x(4e 2t + ^ - 1) 3/2 . We have 


U u 



2 e A sin 





(28) 


where 

4+ = -4- = |e 2t + - 1. (29) 

Equations (19), (23), (26) and (28) give the expression of in regions l)-4), respectively. Note that the term 
oc w 2 in (29) must be discarded in the limit t —> oo, since it gives an artificially large contribution due to the use of 
Pmod{x) instead of P{x) in our analysis above. 


IV. TRANSITION TO QUASI-CLASSICALITY 

Equation (7) describes an oscillator with a time-dependent frequency. Its form is analogous to the form of the 
equation used, for example, to study the evolution of a test massless scalar quantum field in an expanding Universe, 
see e.g. [11], [10] [33], although potential V c {rj) used in the cosmological setting has another time dependence[34], 77 is 
the so-called ’conformal time’. In a class of cosmological models, notably in the important case of inflationary models 
Vc{if) tends to zero at small values of 77 , while at sufficiently large 77 its absolute value can be much larger than the 
square of the wavenumber, fc 2 , characterising a particular field mode. In this case the in -vacuum state specifying initial 
conditions for the field evolution is typically defined in a way analogous to what is used in this paper. Modes evolving 
in the regime k c < \/\V c (if)\ during a considerable period of time experience a copious growth of occupation numbers 
with respect to a natural out-vacuum state. When inflationary models are considered, such a regime corresponds to 
physical wavelengths of the modes being larger than the cosmological horizon scale, and, accordingly, k c C 1/f). In 
the framework of these models it has been shown (see e.g. [11], [13]) that a large scale part of the field consisting of 
the modes with sufficiently small wavenumbers, which have been evolved outside the cosmological horizon for a long 
time behave as a classical random field with zero expectation value. We suppose similar behaviour in our case for 
the part of the wavefunction consisting of modes with u> < at large values of our time variable t. To probe the 

classical behaviour we are going to calculate explicit values of the Bogolyubov coefficients given by expressions (17) 
and, following [13], the time dependence of the Wigner function corresponding to our quantum state, which has a 
special form for a system evolving in a quasi-classical regime. 


A. Calculation of the amount of produced universes 


The amount of universes produced at late ’times’ t —> 00 per unit of volume in k space [35] is given by the square 
of an absolute value of Bogolyubov coefficient /?, 


m 2 = (ft(io)tw + - 1/2, 


(30) 


where we use ( 8 ) and (17). Note that the physical meaning of multiple production of universes was discussed in e.g. 
[ 6 ]. As discussed in this paper |/3 | 2 may be considered as being proportional to a probability to find a universe with 




FIG. 2: The square of absolute value of the Bogolyubov coefficient f) w is shown as a function of ratio uj/uj C rit , for A = 3 * 10 3 
Solid and dashed curves represent numerical result and expression (31), respectively. 


given parameters, in our case characterised by different values of k. In quantum cosmology the picture of multiple 
universes itself has a direct physical meaning either when their wavefunctions can interfere (e.g. [6]) or through 

non-linear interaction terms added to the quantized WdW equation, see e.g. [22]. 

In order to obtain |/3| 2 explicitly we assume that t = to > oo, thus calculating occupation numbers with respect to 

our out-vacuum state. In this limit we have fl(f 0 ) ~ y^e 3 * and (j>+ ~ §e 2t . Substituting these expressions into (28) 
and the resulting one into (30) we obtain 


I/M 2 


1 

8 sinh(^) 



(31) 


[36]. Equation (31) tells that \/3 ul \ 2 oc when oj -C u> C rit ■ This result was obtained in [7] by qualitative means. On 
the other hand, as it is evident from (31) \(3 U \ 2 strongly decreases with increase of w. It is formally equal to zero at 
w ~ However, this value is slightly larger than co cr it and we derived (31) assuming that ui <C w cr jt- Therefore, 
the latter effect may be an artifact of our approximations. Fig. 2 shows |/M 2 given by equation (31) as well as the 
result of calculation of |/3 W | 2 by numerical means. One can see that there is a very good agreement between analytic 
and numerical results at small values of w. When to ~ ui cr it numerical and analytic curves differ by several orders of 
magnitude, but have qualitatively similar behaviour. In both cases there is a drastic decrease of |/M 2 towards u> cr it- 
When w > uj cr it the numerical curve shows that |/3^| 2 is small. 


B. Wigner function 

As was discussed in e.g. [14], [13] in a quasi-classical regime the Wigner function has a special form of a distribution 
over a generalised coordinate multiplied by a sharply peaked distribution centred at a linear combination of generalised 
momentum and coordinate of a classical trajectory. Let us calculate it for our model. In this Section, we assume, for 
simplicity, that our quantum field H/t is quantized in a box, since this assumption doesn’t influence our conclusions. 
This has an advantage that we shall deal with ordinary functions rather than functionals when treating wavefunctions 
corresponding to the quantum state of our system[37]. In this case the Kronecker deltas are implied in (2) and (8) 
instead of delta functions, and summation instead of integration in the expressions containing integrals over k is to 
be used. The normalisation condition also changes but we are nonetheless going to use equation (8) since it doesn’t 
affect our results. 

After this assumption is made, the expression (13) tells that our model is reduced to a discrete set of oscillators 
with time dependent frequency and different values of k and a. Let us consider one of them with a particular value 
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of u>. The Wigner function corresponding to this oscillator has the standard form 

1 r°° 

W(p,q) = - dy^* s {q + y)^ s {q-y)e pv , (32) 

^ J — oo 


where c-numbers p and q are conjugated momentum and coordinate, 'hg(g) is the wave function in the Sclirodinger 
representation, and we omit indices fc, a and u> in this Section. 

To find the evolution of 4 's(q) let us note that in the Schrodinger representation the operators p and q entering 
(13) have the standard form p = — and q = q. On the other hand 4 , g(g) corresponds to the vacuum state, and, 
therefore, c^s(q) = 0. Setting O(to) = w in equations (14)-(16) we express c in terms of p and q to get 


— +Dg ) 4' s (g) = 0, 


where 


D = u> 



(33) 


(34) 


Note that the same equation was obtained in [13] in the quantum field context. Taking into account that when 
t —> —oo 4 ’s(q) tends to that of the vacuum state of an oscillator with frequency u> we have 

’Mg)=(|^) K = D + D*, (35) 

and, substituting (35) into (32) we obtain 


W(j>,q) = - exp 
7r 


Kg 2 (2 p + Rg) 2 \ ' 

2 2K ) y 


R = i(D* — D). 


(36) 


The Wigner function has an important symmetry with respect to interchange of the canonical variables. Namely, it 
can be written in another form 


W (p, q) = - exp 
7r 


KnP 2 (2 q - R n p) 2 \ 

2 2K n ) \ ’ 


K„ 


1 

Jr 


l 

D' 




(37) 


This property follows from the observation that when the momentum rather than coordinate representation for the 
operators p and q is used we have a similar equation for wavefunction differing from (33) by a new coefficient D n = 1/D. 
Also, in the momentum representation there should be the opposite sign of the argument in the exponent in (32). 
Note that since the Wigner function is positive definite in our case we treat it as determining a probability distribution 
in phase space. 

When either K — > 0 or I\ n 0 distributions (36), (37) are sharply peaked around 2p + Rq = 0 and 2 q — R n p, 
respectively. Let us show that these relations hold “on average “ for a classical solution of our equations of motion 
in the limit f > 1 and calculate coefficients K , K ni R and R n explicitly, in the same limit. To do so we simplify the 
expressions for a and (3 in equation (17) using the same approximations as in Section IV A, but now setting Q(to) = oj 
there. We obtain 

a = (2 tt)^~C uj ^(w 1 / 2 fl _1 ^ 2 (t)ai + w -1 / 2 0 1 / 2 (f)a2) sin</> + *(w 1 / 2 0~ 1 / 2 (t)a2 + <W 1,/2 fl 1 ^ 2 (t)ai) cos , (38) 

and 

j3 = (2 tt ) ^ C* ^(w 1 / 2 fl“ 1 / 2 (f)ai — w _1//2 0 1 / 2 (f)a2) sin<^ + «(w 1 / 2 n _1 / 2 (f)ai — w^ 1 / 2 0 1 / 2 (t)a2) cos <j?J , (39) 


where, by definition, f2(f) = y ^ e 3t is the asymptotic value of (15) with to being substituted by t , in the limit t —$■ oo, 

Oi = 2e^, 02 = sinh(^)e - ^, (j> = |a 72 -e 3t (l — |-e _2t ) 3 / 2 + f is the phase of sine and cosine entering (28). Note that 
4> « f2(f) and we use equation (29) for <j> + setting ui = 0 there due to the reason explained above, but taking into 
account the curvature term in the arguments of sine and cosine. 

Substituting (38) and (39) in (34)-(37) we have 


K = 20 


aio 2 


a? sin 2 ( 


20 


a 2 1 
ai sin 2 


R = O 


(a 2 — a\) sin 20 
a? sin 2 6 + ai, cos 2 


—20 cot 0, 


(40) 
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and 


K n = 2 n 


-l 


aia 2 


' + aj cos 2 


2ST 


_1 012 


a i cos z 


R n = Q 


-l 


(af — a 2 ) sin 20 
2 sin 2 0 + a 2 cos 2 


2fl 1 tan 0, 


(41) 


where Q = fl(f) from now on, we take into account that when to <C co C rit we have a 2 <C ai to get the approximate 
expressions. Note that these are valid only when the phase 0 is not very close to <J)j = 7 vj and <j>j = 7r/2 + 7 vj (j is an 
integer), respectively, for equations (40) and (41). Hereafter, we assume that the value of time variable t is such that 
this condition is satisfied. The opposite case should be treated separately. 

When the classical motion is considered p = q. Using this fact and the approximate expressions for R and R n in (40) 
and (41) we easily see that both combinations 2 p + Rq and 2 q — R n p are equal to zero provided that q is proportional 
to sin 0, which is an approximate solution of the classical equations of motion at large times. Thus, in both cases 
of small K and K n distributions (36) and (37) can be treated as describing bundles of classical trajectories with a 
random Gaussian distribution of the coordinate and momentum, respectively. It is interesting to point out that these 
trajectories have the same phase 0. Thus, in the quasi-classical limit only one of two linearly independent solutions of 
the equations of motion is present. This is analogous to the test quantum field problem, where the so-called ’growing’ 
mode is singled out in the same regime. Let us stress again that this picture is not correct when the degenerate set 
of phase 0j is considered. In the latter case the behaviour of the system is no longer quasi-classical. 

i —si 

It is instructive to introduce a ’nascent delta function’ S € (x) = ^= e~~sr having the property that d e _>. o(x) 5(x) 
and rewrite (36) and (37) in an equivalent form as 


W(p,q) 






(42) 


Equation (42) tells that when K —> 0 or K n —> 0 the probability distribution over q or p has a dispersion much 
larger than 1, while the probability distribution over the other coordinate shrinks. This is a characteristics of a 
strongly squeezed vacuum state. It is well known in this case the coordinate having a large dispersion can be treated 
as a classical random quantity with a Gaussian distribution [14].Thus, we have two possible cases of quasi-classical 
evolution of one particular mode referred to hereafter as case one and case two. In the former case corresponding 
to modes with relatively small u> the coordinate can be treated as a classical random variable, while in the latter 
case valid for modes with relatively large u>, which are, however, smaller than uj cr it it is the momentum, which is 
quasi-classical. 

From the condition -U = i J dq i-T = 1 let us find a typical frequency, w s , which separates modes evolving in the 
two different quasi-classical regimes. When ui <C uj s we have case one, while the opposite case two is realised when 
w s < w < oj cr u . Note that when ui ~ io s the system behaviour is essentially quantum. We have 


w s = 


2 

7T 


| - 3t - \ In 
A 2 



(43) 


where we assume, for simplicity, that > 1. Neglecting the logarithmic correction we see that the case one is 
present in the system only when 


t tcrit — 


3A’ 


(44) 


and, accordingly, the scale factor a < a cr n = esx. 


V. QUASICLASSICAL STOCHASTIC WAVEFUNCTION OF THE UNIVERSE AND ITS 
DESCRIPTION THROUGH A DENSITY MATRIX 

A. Quasiclassical wavefunction 

As we discussed in the previous Section when 1 <f < t cr it the large scale part of the T consisting of modes with 
u> <C ui a , Tqc, may be considered as a quasi-classical one. It can be represented using an approach analogous to the 
so-called coarse graining procedure frequently used for quantum fields evolving in an inflationary Universe [11] as an 
integral over all modes with frequencies smaller than lol = where a constant £ l is assumed to be small. It is 

easy to see that wl is given by the same expression as (43), but with the argument of the logarithm being divided by 
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where all quantities are assumed to be c-numbers and the mode function is given by (28). Complex random numbers 
a k and b k have Gaussian distributions and must be normalised in such a way that resulting statistical averages of 
different products of 'I ' qc and its complex conjugate coincide with vacuum averages of products of 4' and 'LL 

< a k a* k , >= S n (k - k), < b k b* k , >= S n (k - k), (46) 

where < ... > denotes a statistical average from now on and all other correlators are equal to zero. 

Let us calculate the averaged square of the absolute value of d ' qci V =< > giving an average probability 

to find a universe. In the beginning we formally assume that lul ~ u> cr u■ Clearly, in this limit V* = V(ui s = u> cr it) = 
(0| 'L'l't |0). We have 


V * = 


-n/2 


e 3t I\ sin 2 </> + e * I 2 cos 2 cj^j , 


where 


h = 


dto 


0 sinh(^) 


2 n T(n/2) 

-1 2" — 1 f 

^7T=2——(n-l)!C(n), h = 

~T) 71 Jo 


dkJUJ r 


^sinhf— 
V 2 


n— 1 


(47) 


(48) 


where we use the same approximations as in the previous Section, (]{n ) is the Riemann zeta function and the integration 
limit in the first integral is formally extended from co cr it to infinity since the integrand there decreases exponentially 
with oj. Remembering that oj cr it = we see that the second term in (47) is much smaller than the first and can 
be discarded. Note that the integrals in (47) are logarithmic diverging when n = 1. We leave this special case for a 
future paper and assume from now on that we have more than one massless field in our model, n >2. Neglecting the 
second term we can rewrite (47) in the form 


V* = C n e A 3t sin 2 4>, C n 


(2"~1) (n-l)!C(n) /3\ 1/2 

2 n-3 7r |n r(f) \AJ 


(49) 


Remarkably, the expression (49) is the same as the probability distribution obtained with help of the Hartle-Hawking 
wavefunction [4], for a closed Universe with Lambda term, see e.g. [2], his equation (8.63) with V(<p) = 

Now we are going to estimate a correction accounting for lul being smaller than LO C rit, but much larger than unity. 
In this case the upper limit of integration in R should be ujl- We have I\ = f^ L duj -^ = fo° ^ sinh^) ~ 
where A I has the same integrand, but with integration being performed from oj s to co cr i t . Approximating sinh(a;) 
there as \e x we obtain 


A / 


2 n+1 


ir n 




-— 3 1 + In 

A 



(50) 


and 


7> = 7\(1-A I/h) 


(51) 


Equations (50) and (51) tell that when t <C t cr u the correction determined by the dependency of lol on time is 
exponentially small and the averaged quasi-classical wavefunction in our Universe is that of Hartle and Hawking. 

Neglecting the contribution to (47) proportional to I 2 corresponds to discarding the term proportional to cos (f) in 
(28). In this case equation (45) can be written as a product of a time dependent regular factor and a stochastic 
function over the spatial variables 


qc = 




d n k (a 


Ak-x _ 

e ak 


C*t 


- lk ' x b* k ) 


(52) 


In the end let us briefly discuss the opposite case t t cr i t . In this case almost all modes except those with very 
low frequencies are in the regime of high frequency quasi-classical behaviour called ’case two’ in the previous Section. 
In this regime, the modes momenta are quasi-classical, and, accordingly, the time derivative of the wavefunction can 
be treated as a random classical variable, while the wavefunction itself is essentially quantum. Differentiating (52) 
over the time we obtain an explicit expression for the time derivative 

2 /AX 1 / 4 f 

^ qc = — l-\ cos 0 J d n k (C„e ikx a k + C*e~ ikx b* k ) (53) 

Note that equation (49) tells that a probability to find a universe with the scale factor larger than critical is expo¬ 
nentially damped on average. On the other hand from equation (53) it follows that the average of 4/ gc always grows 
with time. 
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B. Density matrix 



FIG. 3: We show the dependence of the ratio p(x,x')/F 2 as a function of the relative distance r = \x — x'\. Solid, dashed, 
dotted and dot-dashed curves are calculated for n = 2, 3, 4 and 5, respectively. 


Let us show that the stochastic wavefunction introduced above allows for an alternative description in terms of a 
density matrix with c-number valued matrix elements. For simplicity we assume below that we are in the regime when 
the wavefunction is approximately classical and, accordingly, 1 <C a -C a cr it . Also, we use sometimes hereafter the bra 
and ket notation, for example representing (52) as Hi qc = ( x\ |’F). Let us stress that we are going to use both position 
representation associated with the field amplitudes, as well as the adjoint momentum representation, which are 
clearly unrelated to the position and momentum representations considered in Section V A, where the dynamics of a 
mode amplitude of the wavefunction is treated. 

In order to calculate explicitly the density matrix p(x,x') corresponding to (52) it is convenient to use temporarily 
the momentum representation taking into account that (fc| |x) = ^ 27r | n / 2 e~ lkx . Making the Fourier transform of (52) 
we get an expression similar to (10) with the exception that now the wavefunction in the momentum representation 
is considered as being c-number valued: 

^ k = (k\\^} = (2n) n / 2 F(a)(C u a k + C:b. k ), F(a) = -j= (jj e^sin^, (54) 

where we remind that a k and b k are classical complex Gaussian random numbers with correlation properties given by 
(46). It is convenient to represent these in terms of real random quantities as 

a, k = i(au + <23 + i(a4, - a 2 )), = i(ai - «3 + i(a 4 + 0 : 2 )), < a t a k >= 6 ik S n (k - k'), (55) 

introduce real and imaginary parts of ^ k and as R = i?e(4'/ c ), / = /m('Ffc), A u = Re{C u ) and B u = Im(Cuj), 
and represent (54) in the form 

R = {2^) n / 2 F{a){A u a l + B u a 2 ), I = (2 K) n/2 F{a){A u a A + B u a 3 ). (56) 

Now, since both R and / are sums of two uncorrelated random Gaussian numbers the general theorem tells that their 
distribution functions V(R) and V(I) are also Gaussian, with the square of dispersion a 2 =< RR >=< II >, while 
their joint distribution function V(R,I) = V(R)V(I). Explicitly, we have 

7W) = 2^ 6XP (~^f)’ = (2n) n F 2 \C u \ 2 . (57) 

The expression (57) gives a probability to find ty k with particular values of R and /. Using this fact the density 
matrix in the momentum representation, p(k,k'), can be obtained by the usual rule p(k,k') = f dRdl’i' k A? k /V(R, /), 
where we can set k = k' under the integral taking into account that p(k , k') is clearly diagonal in this representation. 
In this way we obtain 


p(k, k’) = 2(2Tr) n F 2 \C u \ 2 6 n (k - k'), \C U \ 2 



(58) 
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where we use (19). The expression (58) tells that mainly the modes of with small values of ui give contributions 
to the mixed state defined by p(k , k'), while the contribution of modes with uj 1 is exponentially damped, see also 
the end of this Section. 

The density matrix in the position representation, p(x,x') = / 2 n) n / d n kd n k'e l ^ kx ~ k x ^ p{k,k') is easily obtained 
from (58) taking into account that we can extend the upper limits of integrals over k and k' to infinity, since the 
contribution of modes with w > w s to them is negligible provided that co s 1. We have 

2 n . / ttuj \ 

p(x,x)= dss*Jn=2(rs) sinh - (—J , (59) 

where r = \x — x'\, J^(z) is the Bessel function and we use the standard result that a multidimensional Fourier 
transform of a function depending only on the modulus of the momentum vector takes the form of one-dimensional 
Hankel transform. Note that, since the Bessel function in (59) has a half-integer index, it can be expressed through 
elementary functions. 

That the density matrix is a function only of the relative distance in the field space is clearly due to the homogeneity 
of our model with respect to the field coordinates, which, in its turn, is valid only for massless fields. In general, 
for the potentials V depending on the field coordinates the homogeneity is broken and the density matrix must be 
determined by the field coordinates themselves. The dependence of p{x,x') on r is shown in Fig. 3 for several values 
of n. 

The diagonal elements of the density matrix p(x , x) give probabilities to find fields with particular values of x. They 
can be easily obtained from (50) using the decomposition of Bessel function near the origin of its argument and the 
first expression in (48). In this way we obtain p(x,x) = V*, where V * is given by equation (49) as it should be. 



FIG. 4: The dependence of the expectation values given by equation (61) is shown as a function of n. Solid, dashed, dotted 
and dot-dashed curves are calculated for k = 2, 3, 4 and 5, respectively. 


It is well known that the Fourier transform of (59) determines a probability distribution of particles (universes, 
in our case) over k. and, accordingly, over the field velocities <pi, see e.g. [23]. Then equation (58) tells that it is 
proportional to |C(u;)| 2 and we employ this result to obtain a normalised probability distribution over w, ’P(w), as 


v M = 


UJ 


L |CV)| 2 


1 sinh -1 '' 7ruM 


(T) 


/ 0 °° dujuj n 1 |C(w)| 2 f^° duw" -1 sinh 1( 


( 2 ^) 2(2" - l)(n - l)!C(n) 


sinh 


-(¥) 


(60) 


where we use (48) and (58). Let us remind that in the classical limit w is proportional to the absolute value of 
the ’total’ field velocity, ui oc v to t = \/Xa=i Thus, equation (60) may be used to find probabilities of universes 
with different Vtot- This, in its turn, can be employed to determine natural initial conditions for classical evolution 
of the Universe (or its sufficiently homogeneous parts, see Discussion below). It is instructive to calculate different 
expectation values of powers of oj 


1 (n + fc - 1)! (2 n+k - 1) C (n + k ) 

7 r k (n — 1)! (2™ — 1) C( n ) 


( 61 ) 


where we use again (48). The dependence of < u k > on n is shown in Fig. 4 for k=2..5, respectively. 
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VI. CONCLUSIONS AND DISCUSSION 

We show, in the framework of a simple thirdly quantized minisuperspace model of a closed FRW universe with 
a small Lambda term and n massless scalar fields that its wavefunction operator has a simple interpretation in the 
limit of large scale factors a 1 provided that a natural m-vacuum state is specified for the system. Namely, when 
1 € « < e® the wavefunction operator may be approximately treated as a classical random field with its averaged 
value being proportional to the Hartle Hawking wavefunction. When a it is the derivative of the wavefunction, 

which has the property of a classical random field [38]. 

The physical explanation of this result is the same as in the well developed theory of creation of a test field or 
density inhomogenities/gravitational waves in a inflationary Universe. Both models have a copious production of 
excitations with respect to a suitably chosen out-vacuum state. In both models the in-vacuum in the Sclirodinger 
picture evolves to a strongly squeezed state, while the Wigner function of a particular mode has a special form at large 
values of the natural time variable (or, the scale factor) being proportional to a product of a Gaussian distribution 
of one of the canonical variables with large dispersion and a distribution strongly peaked about a relation involving 
canonical variables, which is valid on solutions to classical equations of motion. 

Thus, the natural in-vacuum for the quantum WdW model provides a well defined classical although stochastic 
wavefunction at sufficiently large values of the scale factor. The fact that its averaged value is proportional to the 
Hartle Hawking result needs a further investigation. Perhaps, some progress could be made in the path integral 
formulation of our theory, although of course the classical action entering the Hartle Hawking formalism should be 
substituted by the action induced by Lagrangian (4) in our case. 

We also show that the stochastic wavefunction formalism is equivalent to the presence of a density matrix describing 
a mixed state. Similar to the wavefunction approach, the matrix elements can be treated as c-numbers at sufficiently 
large values of scale factors. Thus, in the framework of our model a mixed state defined over classical solutions of 
WdW equations emerges in this asymptotic limit. Diagonal elements of the density matrix giving a probability to 
find some particular values of fields and scale factor are again proportional to the corresponding expression following 
from the Hartle and Hawking formalism, thus the density matrix has a trivial form in the position representation. 
However, in the momentum representation the dependence is non-trivial, it gives a probability to find universes with 
different values of field velocities. We calculate this as well as associated expectation values in Section VB. It is 
important to stress that this result is different from what follows from the Hartle-Hawking wavefunction, since the 
latter is uniform in the position representation for the model with massless fields, and, therefore, predicts that field 
velocities are strictly zero in this case. 

Clearly, that the density matrix gives a non-trivial distribution in the momentum representation is due to the fact 
that in our models the field velocities have a direct physical meaning, while the field amplitudes are defined up to 
arbitrary constant values. In principal, the distributions of this kind can be used to specify the most natural initial 
conditions for classical evolution of the Universe. 

Note, however, that program is hampered by the usual problems with interpretation of wavefunction and meaning 
of measurements in quantum cosmology. Indeed, for example, it would be difficult to probe such distributions without 
a consideration of an ensemble of universes with observers belonging to different universes able to communicate with 
each other. Of course, it is difficult to achieve even when we take into account that the universes could ’interact’ 
through quantum interference. The difficulty can be alleviated either in models, where there is a non-linear interaction 
among the universes or in a modification of the model, where it is assumed that it describes sufficiently large locally 
homogeneous parts of positive curvature of a generally inhomogeneous Universe. In the latter case measurements can 
be performed in all parts and compared with each other after, say, the Lambda term decays due to some reason and 
these parts come into causal contact with each other. Similar approaches were recently discussed in e.g. [25] in the 
framework of loop quantum gravity and [26] in a thirdly quantized model. 

It is interesting to point out that when the Universe is assumed to be slightly inhomogeneous, there is another 
mechanism of emergence of a mixed state by averaging out inhomogeneous degrees of freedom, see e.g. [27] or [28]. 
Also, let us note that the density matrix approach has been considered in the third quantization formalism, although 
in a sense quite different from what is discussed in this Paper, see e.g. [29]. 

We suspect that a similar behaviour of wavefunction (or, emergence of the ’classical’ density matrix) at large values 
of a is present in all models, where a large production of excitations is observed. 

Note that, strictly speaking, our results are valid for the number of massless fields n > 2. When n = 1 the integrals 
I\ and I 2 in equation (48) experience a logarithmic divergency at small values of to. This is also similar to the 
well known logarithmic infrared divergency of test field in inflationary models, see e.g. [10] and references therein. 
Although this case definitely needs an additional study, our formalism may be used to define various conditional 
probabilities. 

Finally, it would be very important to generalise our approach to the much more realistic case of a massive scalar 
field in a FRW universe. In this case the separation of variables in the WdW equation is absent. However, this 
equation can be tackled numerically and our approach could be used to interpret numerical results. It is expected 
that the massive field case could have a non-trivial dependence of the density matrix on the field coordinates. 



15 


Acknowledgments 

We are grateful to D. A. Kompaneets, A. A. Starobinsky and V. N. Strokov for valuable discussions, D. M. Gangardt 
and V. N. Lukash for useful remarks and Roma Basko for stimulating questions. This study was supported in part 
by the grant of the President of the Russian Federation for Support of Leading Scientific Schools NSh-4235.2014.2. 


[1] B. S. DeWitt, Phys. Rev. 160, 1113 (1967). 

[2] C. Kiefer, ’Quantum Gravity’, second edition, Oxford University Press (2007). 

[3] A. Vilenkin, Phys. Rev. D 30, 509 (1984). 

[4] J. B. Hartle, S. W. Hawking, Phys. Rev. 28, 2960 (1983). 

[5] T. Banks, Nucl. Phys. B 309, 493 (1988). 

[6] A. Hosoya, M. Morikawa, Phys. Rev. D 39, 1123 (1989). 

[7] V. A. Rubakov, Physics Letters B 214, 503 (1988). 

[8] N. D. Birrell, P. C. W. Davies, “Quantum fields in curved space”, Chicago University Press, (1982). 

[9] A. D. Linde, Physics Letters B 175, 395 (1986). 

[10] A. D. Linde, “Particle Physics and Inflationary Cosmology”, Harwood Academic Publishers, (1990). 

[11] A. A. Starobinsky, Lecture Notes in Physics 246, 107 (1986). 

[12] C. Kiefer, D. Polarski, Annalen der Physik 7, 137 (1998). 

[13] D. Polarski, A. A. Starobinsky, Classical and Quantum Gravity 13, 377 (1996). 

[14] L. P. Grishchuk, Yu. P. Sidorov, Phys. Rev. D 42, 3413 (1990). 

[15] J. J. Halliwell, Phys. Rev. D 39, 2912 (1989). 

[16] Yu. V. Sidorov, Europhys. Letters 10, 407 (1989). 

[17] S. P. Kim, Nuclear Physics B 246, 68, (2014). 

[18] S. Abe, Phys. Rev. D 47, 718 (1993). 

[19] T. M. Dunster, SIAM J. Math. Anal. 21, 995 (1990). 

[20] L. P. Grishchuk, JETP 40, 409 (1975). 

[21] V. N. Lukash, JETP 52, 807 (1980). 

[22] M. McGuigan, Phys. Rev. D 38, 3031 (1988). 

[23] L. P. Pitaevskii, E. M. Lishitz, “Statistical Physics, Part 2: Theory of Condensed State”, Vol. 9, 1st ed., Butterworth- 
Heinemann (1980). 

[24] R. Colistete, R., J. C. Fabris, N. Pinto-Neto, Phys. Rev. D 57, 4707 (1998). 

[25] M. Bojowald, A. L. Chinchilli, D. Simpson, C. C. Dantas, M. Jaffe, Phys. Rev. D 86, 124027, (2012). 

[26] G. Calgani, S. Gielen, D. Oriti, Classical and Quantum Gravity 29, 105005 (2012). 

[27] T. Padmanabhan, Phys. Rev. D 39, 2924 (1989). 

[28] M. Morikawa, Phys. Rev. D 40, 4023 (1989). 

[29] M. A. Castagnino, A. Gangni, F. D. Mazzitelli, I. I. Tkachev, Classical and Quantum Gravity 10, 2495 (1993). 

[30] A care should be taken with this terminology. Rather, as in the general case of creation of particles by an external field, see 
e.g. [8], we are talking about amplification of vacuum fluctuations of the quantized wavefunction by evolving gravitational 
field. 

[31] The fact that the vacuum state is strongly squeezed in a similar model was also note by [16]. 

[32] Note that often the logarithm of scale factor is denoted as a, see e.g. [6], [18]. However, we prefer to use t for this variable 
to stress its time-like character. 

[33] Note that equations describing the evolution of scalar and tensor cosmological perturbations also have the same structure, 
see e.g. [21] and [20]. 

[34] Note that in the cosmological context the potential is often defined with the opposite sign. 

[35] In the classically allowed region the phase S plays the role of classical action characterising a particular classical trajectory of 
our model. Using this fact, one can show that every component of vector k is proportional to the derivative of corresponding 
field <j>i over the physical time. 

[36] In the course of the calculation we use the well known relation T(1 + i*)r(l — ix) = ■ 

[37] Let us stress the difference between a wavefunction corresponding to a quantum state of our system and the wavefunction 
of the Universe. The latter plays the role of a generalised coordinate in the formalism of third quantization. 

[38] Let us note that in certain minisuperspace models considering classical WdW wave functions they also demonstrates 
quasiclassical (in this case WKBJ-like) behaviour only at sufficiently small values of scale factors, see e.g. [24]. 



